Data prep

Cleaning and formatting

We start by preparing the data from the outputs of scripts 2.6, 7.1 and 8.1. The results table shown as output summarised all simulation efforts for the thesis.

# Land use data 
trans_df <- read_csv("../outputs/final/land_use_trans_df.csv")
vals_df <- read_csv("../outputs/final/land_use_vals_df.csv")

# Results dir + mun and mrc shapefile
mun <- st_read("../data/mun/munic_SHP_clean.shp", quiet = TRUE)
mrc <- st_read("../data/raw/vector/mrc_SHP/mrc_polygone.shp",  quiet = TRUE)
mrc.mont <- mrc %>% filter(MRS_NM_REG=="Montérégie") 
mrc.mont.reproj <- st_transform(mrc.mont, raster::crs(mun))

# Classes
classes <- read_csv("../config/rcl_tables/land_use/recode_table.csv")
forest_classes <- read_csv("../config/rcl_tables/land_use/recode_table_forest.csv") %>% 
  rename(new_code = "ID", new_class = "Name")

classes_unique <- unique(classes[, c("new_code", "new_class")]) %>% 
  bind_rows(forest_classes)

# results key
sce_code_vec <- 
  as.vector(t(outer(c("BAU-", "R-", "CorPr-", "R-CorPr-", "R(T)-CorPr-"), 
                    c("Hist", "Base", "RCP8"), paste0)))
sce_code_vec_run <- 
  rep(c("BAU", "R", "CorPr", "R-CorPr", "R(T)-CorPr"), each=3)

results_clean <- read_csv("../outputs/final/stsim_run_results.csv") %>% 
  dplyr::select(scenarioId, name) %>% 
  mutate(name =  gsub(.$name, pattern = " \\(.+\\)", replacement = "")) %>% 
  mutate(sce =  paste0("sce_", .$scenarioId)) %>% 
  mutate(chapter = c("none", rep("both", 3), rep("chap_1", 3), rep("chap_2", 9))) %>% 
  mutate(splitted = str_split(name, " \\| ")) %>% 
  mutate(climate = unlist(map(splitted, ~unlist(.x[2])))) %>% 
  mutate(run = unlist(map(splitted, ~unlist(.x[1])))) %>% 
  dplyr::select(-splitted) %>% 
  replace_na(list(climate = "none")) %>% 
  mutate(code = c("Control", sce_code_vec)) %>% 
  mutate(code_run = c("control", sce_code_vec_run))

# Final extracted datasets
df_final <- read_csv("../outputs/final/final_df_current_density_part1.csv") %>%
  bind_rows(read_csv("../outputs/final/final_df_current_density_part2.csv")) %>% 
  mutate(timestep = (timestep*10)+1990, source = "model", zone = as.character(zone))
df_final_origin <- read_csv("../outputs/final/final_df_origin_current_density.csv") %>%
  mutate(timestep = timestep*10+1990, source = "model")

# Stratum key
stratum_key <- read_csv("../config/stsim/SecondaryStratum.csv") %>%
  mutate(ID = as.factor(ID)) %>%
  rename(zone=ID, MUS_NM_MUN=Name) %>%
  mutate(zone = as.character(zone)) %>% 
  left_join(df_final, by="zone") %>%
  filter(MUS_NM_MUN!="Not Monteregie")

# Sum all municipalities
df_summarised <- df_final %>%
  group_by(sce, timestep, species, iteration) %>% 
  summarise(sum_cur = sum(current)) %>% ungroup %>%
  mutate(source = "model")
df_origin_summarised <- df_final_origin %>%
  group_by(timestep, species) %>% 
  summarise(sum_cur = sum(mean)) %>% ungroup %>% 
  mutate(sce = "sce_0", source = "observation")

# Joined dataset
joined <- full_join(df_summarised, df_origin_summarised, 
                    by=c("sce", "source", "species", "timestep", "sum_cur")) %>% 
  left_join(results_clean, by = "sce")  %>% 
  replace_na(list(climate = "none", run = "historic run")) %>% 
  # Make factors
  mutate(sce = as.factor(sce), run = as.factor(run), sum_cur = 10*(sum_cur),
         climate = factor(climate, levels = c("none", "historic",
                                              "baseline", "RCP 8.5")),
         run = factor(run, levels = c("historic run", "BAU run", 
                                      "BAU run + corrs protection", "BAU run + ref", 
                                      "BAU run + corrs protection + ref",
                                      "BAU run + corrs protection + ref TARGETED")))

# Histogram data 
hist_original <- read_csv("../outputs/final/final_values_output_original.csv")
hist_true <- read_csv("../outputs/final/final_values_output_TRUE.csv") %>% 
  mutate(sce = "TRUE")
histograms <- read_csv("../outputs/final/final_values_output.csv") %>% 
  bind_rows(hist_original) %>% 
  bind_rows(hist_true) %>% 
  left_join(results_clean, by="sce") 

# SURF Data
surf <- read_csv("../surf/surf_output.csv") %>% 
  mutate(timestep = timestep*10+1990) %>% 
  left_join(results_clean, by=c("scenario"="scenarioId")) %>% 
  replace_na(list(climate = "none", run = "historic run")) %>% 
  # Make factors
  mutate(sce = as.factor(sce), run = as.factor(run),
         climate = factor(climate, levels = c("none", "historic",
                                              "baseline", "RCP 8.5")),
         run = factor(run, levels = c("historic run", "BAU run", 
                                      "BAU run + corrs protection", "BAU run + ref", 
                                      "BAU run + corrs protection + ref",
                                      "BAU run + corrs protection + ref TARGETED")))

# Bar plot data
bar_data <- read_csv("../outputs/final/final_bar_plot_data.csv")

results_clean

Figures

Mean current flow change

Historic

the_width = 18 
the_height = 15
the_dpi = 500

fig_1_historic <- joined %>% 
  filter(climate == "none") %>% 
  group_by(scenarioId, species, iteration) %>% 
  group_modify(~mutate(.x , baseline = subset(.x, timestep == 1990)$sum_cur, 
                       the_diff = ((sum_cur-baseline)/baseline)*100 )) %>% 
  ggplot(aes(x=timestep, y=the_diff, col=source)) +
  scale_color_manual(values=c('#d8b365','#5ab4ac'), 
                     labels = c("Model", "Observation")) +
  geom_smooth(aes(group = sce), method = "glm", se=FALSE) +
  scale_linetype_manual(values = c(1,3,2,5))+
  geom_point(aes(group = sce)) +
  facet_wrap(~species, scales = "fixed") +
  labs(y = "Current flow (% Change)",
       x = "Year",
       col = "Source") +
  NULL
fig_1_historic
Connectivity change for species through time, 1990-2010

Connectivity change for species through time, 1990-2010 

ggsave(fig_1_historic,
       filename = "../thesis/figures/connectivity_decrease_x5species_historic.png",
       width = the_width, height = the_height, dpi = the_dpi)

Chap 1

Line graph
colors_sce <- 
  data.frame(sce = c("Hist", "Baseline", "RCP85"), 
             color = c("#8DA0CB", "#FC4E2A", "#800026"), stringsAsFactors = F)

fig_1_static_chap_1 <- joined %>% 
  filter(climate != "none", chapter %in% c("none", "both", "chap_1")) %>% 
  group_by(scenarioId, species, iteration) %>% 
  group_modify(~mutate(.x , baseline = subset(.x, timestep == 2010)$sum_cur, 
                       the_diff = ((sum_cur-baseline)/baseline)*100 )) %>% 
  ggplot(aes(x=timestep, y=the_diff, col=climate)) +
  scale_color_manual(values = colors_sce$color, 
                     labels = c("Historic", "Baseline", "RCP 8.5")) +
  geom_smooth(aes(group = sce, linetype=code_run), alpha=0.5, method="loess") +
  #geom_point(aes(group = sce)) +
  scale_linetype_manual(values=c("solid", "twodash")) +
  facet_wrap(~species, scales = "fixed") +
  geom_hline(yintercept = 0, linetype = "dashed", alpha = 0.5) +
  #geom_vline(xintercept = 2020, linetype = "dashed", color = "darkGreen", alpha = 0.5) +
  labs(y = "Current flow (% Change)",
       x = "Year",
       col = "Climate",
       linetype = "Run") +
  guides(col = guide_legend(ncol = 2, nrow=2), 
         linetype = guide_legend(nrow=1, override.aes = list(colour = 'black'))) +
  theme(legend.position = c(0.85, 0.15))+
  NULL
fig_1_static_chap_1
Connectivity change for species through time, 2010-2100

Connectivity change for species through time, 2010-2100 

ggsave(fig_1_static_chap_1, 
       filename = "../thesis/figures/connectivity_decrease_x5species_chap1.png", 
       width = the_width, height = the_height, dpi = the_dpi)
Radar Graph
# Default size
gridline.label.offset   = 4
legend.text.size = 20
group.line.width = 0.9
grid.label.size = 10
background.circle.colour    = "#ffffff"
group.point.size    = 5
axis.label.size = 8
group.colours = RColorBrewer::brewer.pal(name = "Paired",n=10)[c(2,4,6,8,10)]


radar_data_chap1 <- joined %>% 
  filter(climate != "none", chapter %in% c("none", "both", "chap_1")) %>% 
  group_by(timestep, species, code) %>% 
  summarise(sum_cur = mean(sum_cur)) %>% ungroup %>% 
  pivot_wider(names_from = timestep, values_from = sum_cur) %>% 
  mutate(diff = ((`2100`-`2010`)/`2010`)*100) %>% 
  select(-c(`2010`:`2100`)) %>%  
  pivot_wider(names_from = code, values_from = diff) %>% 
  rename(group = species)

radar_chap1 <-  
  ggradar(radar_data_chap1, centre.y = -20, legend.position = "right",
          grid.min = -20, grid.max = 3, grid.mid = 0, 
          values.radar = c("-20 %", "0 %", ""), 
          group.colours = group.colours,
          gridline.label.offset = gridline.label.offset, 
          legend.text.size = legend.text.size,
          group.line.width = group.line.width, 
          grid.label.size = grid.label.size,
          background.circle.colour = background.circle.colour, 
          group.point.size  = group.point.size, 
          axis.label.size = axis.label.size)
radar_chap1

ggsave("../thesis/figures/radar_ggradar_chap1.png", radar_chap1,
       width = the_width, height = the_height, dpi = the_dpi)

Chap 2

fig_1_static_chap_2 <- joined %>% 
  filter(climate != "none", chapter %in% c("none", "both", "chap_2")) %>% 
  group_by(scenarioId, species, iteration) %>% 
  group_modify(~mutate(.x , baseline = subset(.x, timestep == 2010)$sum_cur, 
                       the_diff = ((sum_cur-baseline)/baseline)*100 )) %>% 
  ggplot(aes(x=timestep, y=the_diff, col=climate)) +
  scale_color_manual(values = colors_sce$color,
                     labels = c("Historic", "Baseline", "RCP 8.5")) +
  geom_smooth(aes(group = sce, linetype=code_run), alpha=0.5, method="loess") +
  scale_linetype_manual(values=c("solid", "dashed", "longdash", "dotted")) +
  facet_wrap(~species, scales = "fixed") +
  geom_hline(yintercept = 0, linetype = "dashed", alpha = 0.5) +
  geom_vline(xintercept = 2020, linetype = "dashed", color = "darkGreen", alpha = 0.5) +
  labs(y = "Current flow (% Change)",
       x = "Year",
       col = "Climate",
       linetype = "Run") +
  guides(col = guide_legend(ncol = 2, nrow=2), 
         linetype = guide_legend(nrow = 2, ncol=2, 
                                 override.aes = list(colour = 'black'))) +
  theme(legend.position = c(0.85, 0.15))+
  NULL
fig_1_static_chap_2
Connectivity change for species through time, 2010-2100

Connectivity change for species through time, 2010-2100 

ggsave(fig_1_static_chap_2, 
       filename = "../thesis/figures/connectivity_decrease_x5species_chap2.png", 
       width = the_width, height = the_height, dpi = the_dpi)
Radar Graph
radar_data_chap2 <- joined %>% 
  filter(climate != "none", chapter %in% c("none", "both","chap_2")) %>% 
  group_by(timestep, species, code) %>% 
  summarise(sum_cur = mean(sum_cur)) %>% ungroup %>% 
  pivot_wider(names_from = timestep, values_from = sum_cur) %>% 
  mutate(diff = ((`2100`-`2010`)/`2010`)*100) %>% 
  select(-c(`2010`:`2100`)) %>%  
  pivot_wider(names_from = code, values_from = diff) %>% 
  rename(group = species)

radar_chap2 <-  
  ggradar(radar_data_chap2, centre.y = -20, legend.position = "right",
          grid.min = -20, grid.max = 3, grid.mid = 0, 
          values.radar = c("-20 %", "0 %", ""), 
          group.colours = group.colours,
          gridline.label.offset = gridline.label.offset, 
          legend.text.size = legend.text.size,
          group.line.width = group.line.width, 
          grid.label.size = grid.label.size,
          background.circle.colour = background.circle.colour, 
          group.point.size  = group.point.size, 
          axis.label.size = axis.label.size-2)
radar_chap2

ggsave("../thesis/figures/radar_ggradar_chap2.png", radar_chap2,
       width = the_width, height = the_height, dpi = the_dpi)

Both

Radar Graph
radar_data_all <- joined %>% 
  filter(climate != "none") %>% 
  group_by(timestep, species, code) %>% 
  summarise(sum_cur = mean(sum_cur)) %>% ungroup %>% 
  pivot_wider(names_from = timestep, values_from = sum_cur) %>% 
  mutate(diff = ((`2100`-`2010`)/`2010`)*100) %>% 
  select(-c(`2010`:`2100`)) %>%  
  pivot_wider(names_from = code, values_from = diff) %>% 
  rename(group = species)

radar_all <-  
  ggradar(radar_data_all, centre.y = -20, legend.position = "right",
          grid.min = -20, grid.max = 3, grid.mid = 0, 
          values.radar = c("-20 %", "0 %", ""), 
          group.colours = group.colours,
          gridline.label.offset = gridline.label.offset, 
          legend.text.size = legend.text.size,
          group.line.width = group.line.width, 
          grid.label.size = grid.label.size,
          background.circle.colour = background.circle.colour, 
          group.point.size  = group.point.size, 
          axis.label.size = axis.label.size-2)
radar_all

ggsave("../thesis/figures/radar_ggradar_both.png", radar_all,
       width = the_width, height = the_height, dpi = the_dpi)

ROC Curves

ROC curves from fitting both models, plotting with patchwork package.

# Roc curves data
urb <- readRDS("../data/temp/fit_rs_outcome_rf_urb_2.RDS")
agex <- readRDS("../data/temp/fit_rs_outcome_rf_agex_2.RDS")

p1 <- bind_rows(urb, .id = "fold") %>% 
  mutate(Fold = factor(fold, levels = as.character(1:10))) %>% 
  group_by(Fold) %>% 
  roc_curve(.pred, truth = outcome_fact) %>% 
  autoplot(add=T) +
  ggtitle(paste("Urbanisation")) +
  annotate(x = 0.75, y = 0.25, geom="label", size = 3,
           label = as.character(paste("Av AUC = 0.938 +/- 0.002"))) +
  theme(legend.position = "none") +
  NULL

p2 <- bind_rows(agex, .id = "fold") %>% 
  mutate(Fold = factor(fold, levels = as.character(1:10))) %>% 
  group_by(Fold) %>% 
  roc_curve(.pred, truth = outcome_fact) %>% 
  autoplot(add=T) +
  ggtitle(paste("Agricultural Expansion")) +
  annotate(x = 0.75, y = 0.25, geom="label", size = 3,
           label = as.character(paste("Av AUC = 0.929 +/- 0.002"))) +
  NULL

full = p1 + p2

full
ggsave("../thesis/figures/double_roc_resample.png", full)

Histograms

Original

bi_colors <- c("#67a9cf", "#ef8a62")

hist_plot_origin <- histograms %>%
  filter(sce %in% c("TRUE", "sce_37")) %>% 
  filter(ts == 2) %>% 
  mutate(ts = ts*10+1990) %>% 
  mutate(ts = as.factor(ts)) %>% 
  ggplot(aes(x=bins, y=n, fill=sce, colour =sce)) + 
  geom_density(stat='identity', show.legend=T, alpha=0.3)+
  scale_fill_manual(values = bi_colors, name="Source", labels=c("Model", "Observation")) +
  scale_color_manual(values = bi_colors, name="Source", labels=c("Model", "Observation")) +
  facet_wrap(~species) +
  labs(x = "Flow intensity distribution (log)", 
       y = "") +
  theme(strip.text.y = element_text(angle=360, size=10, hjust = 0), 
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank(), 
        axis.text.y.left = element_blank(),
        axis.ticks.y.left = element_blank()) +
  NULL
hist_plot_origin

ggsave("../thesis/figures/original_hists.png", hist_plot_origin,
       width = the_width, height = the_height, dpi = the_dpi)

Chap 1

bi_colors <- c("#67a9cf", "#ef8a62")

hist_plot_1 <- histograms %>% 
  mutate(ts = as.factor(ts)) %>% 
  filter(chapter %in% c('both', 'chap_1')) %>%  
  ggplot(aes(x=bins, y=n, group=ts)) + 
  geom_density(stat='identity', show.legend=T, aes(fill=factor(ts), color=factor(ts)), alpha=0.3)+
  facet_grid(code~species) +
  scale_fill_manual(values = bi_colors, 
                    labels = c("2010", "2100")) +
  scale_color_manual(values = bi_colors, 
                     labels = c("2010", "2100")) +
  labs(x = "Flow intensity distribution (log)", 
       y = "")+
  theme(strip.text.y = element_text(angle=360, size=10, hjust = 0), 
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank(), 
        axis.text.y.left = element_blank(),
        axis.ticks.y.left = element_blank(), 
        legend.position = c(0.95, 0.96)) +
  labs(color="Timestep", fill="Timestep") +
  NULL
hist_plot_1

ggsave("../thesis/figures/hist_chap1.png", hist_plot_1, width = 18, height = 15)

Chap 2

hist_plot_2 <- histograms %>% 
  mutate(ts = as.factor(ts)) %>% 
  filter(chapter %in% c('chap_2')) %>%  
  ggplot(aes(x=bins, y=n, group=ts)) + 
  geom_density(stat='identity', aes(fill=factor(ts),  colour=factor(ts)), alpha=0.3)+
  facet_grid(code~species) +
  scale_fill_manual(values = bi_colors, 
                    labels = c("2010", "2100")) +
  scale_color_manual(values = bi_colors, 
                     labels = c("2010", "2100")) +
  labs(x = "Flow intensity distribution (log)", 
       y = "") +
  theme(strip.text.y = element_text(angle=360, size=10, hjust = 0), 
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank(), 
        axis.text.y.left = element_blank(),
        axis.ticks.y.left = element_blank(),
        legend.position = c(0.95, 0.96)) +
  labs(color="Timestep", fill="Timestep") +
  NULL
hist_plot_2

ggsave("../thesis/figures/hist_chap2.png", hist_plot_2, width = 18, height = 20)

SURF Analysis

Chap1

Linear Graph
surf_1 <- surf %>% 
  mutate(scenario = as.factor(scenario)) %>% 
  filter(climate != "none", chapter %in% c("none", "both", "chap_1")) %>% 
  filter(climate != "none") %>% 
  group_by(scenario, species, iter) %>% 
  group_modify(~mutate(.x , baseline = subset(.x, timestep == 2010)$kp_nb, 
                       the_diff = ((kp_nb-baseline)/baseline)*100 )) %>% ungroup %>% 
  ggplot(aes(x = timestep, y = the_diff, color=climate)) +
  scale_linetype_manual(values=c("solid", "twodash")) +
  geom_smooth(aes(group = scenario, linetype=code_run), alpha=0.2, method="loess")+
  #geom_point(aes(group=scenario), alpha= 0.05, size = 0.5)+
  facet_wrap(~species, scales = "fixed")+
  scale_color_manual(values = colors_sce$color,
                     labels = c("Historic", "Baseline", "RCP 8.5")) +
  labs(y = "Change in nb of Keypoints detected (%)",
       x = "Year",
       col = "Climate",
       linetype = "Run") +
  guides(col = guide_legend(ncol = 2, nrow=2), 
         linetype = guide_legend(nrow=1, override.aes = list(colour = 'black'))) +
  NULL
surf_1

ggsave("../thesis/figures/surf_chap1.png", surf_1,
       width = the_width, height = the_height, dpi = the_dpi)
Radar Graph
surf_radar_data_1 <- surf %>% 
  filter(climate != "none", chapter %in% c("none", "both","chap_1")) %>% 
  filter(sce != "sce_0", name != "historic run") %>% 
  group_by(timestep, species, code) %>% 
  summarise(kp_nb = mean(kp_nb)) %>% ungroup %>% 
  pivot_wider(names_from = timestep, values_from = kp_nb) %>% 
  mutate(diff = ((`2100`-`2010`)/`2010`)*100) %>% 
  select(-c(`2010`:`2100`)) %>%  
  pivot_wider(names_from = code, values_from = diff) %>% 
  rename(group = species)

surf_radar_1 <- 
  ggradar(surf_radar_data_1, centre.y = -60, legend.position = "right",
          grid.min = -60, grid.max = 5, grid.mid = 0, 
          values.radar = c("-60 %", "0 %", ""), 
          group.colours = group.colours,
          gridline.label.offset = 20, 
          legend.text.size = legend.text.size,
          group.line.width = group.line.width, 
          grid.label.size = grid.label.size,
          background.circle.colour = background.circle.colour, 
          group.point.size  = group.point.size, 
          axis.label.size = axis.label.size)
surf_radar_1

ggsave("../thesis/figures/surf_radar_chap1.png", surf_radar_1,
       width = the_width, height = the_height, dpi = the_dpi)

Chap 2

Linear Graph
surf_2 <- surf %>% 
  mutate(scenario = as.factor(scenario)) %>% 
  filter(climate != "none", chapter %in% c("none", "both", "chap_2")) %>% 
  filter(climate != "none") %>% 
  group_by(scenario, species, iter) %>% 
  group_modify(~mutate(.x , baseline = subset(.x, timestep == 2010)$kp_nb, 
                       the_diff = ((kp_nb-baseline)/baseline)*100 )) %>% ungroup %>% 
  ggplot(aes(x = timestep, y = the_diff, color=climate)) +
  geom_smooth(aes(group = scenario, linetype=code_run), alpha=0.2, method="loess")+
  #geom_point(aes(group=scenario), alpha= 0.05, size = 0.5)+
  scale_linetype_manual(values=c("solid", "dashed", "longdash", "dotted")) +
  facet_wrap(~species, scales = "fixed")+
  scale_color_manual(values = colors_sce$color,
                     labels = c("Historic", "Baseline", "RCP 8.5")) +
  labs(y = "Change in nb of Keypoints detected (%)",
       x = "Year",
       col = "Climate",
       linetype = "Run") +
  guides(col = guide_legend(ncol = 2, nrow=2), 
         linetype = guide_legend(nrow=2, ncol=2, 
                                 override.aes = list(colour = 'black'))) +
  NULL
surf_2

ggsave("../thesis/figures/surf_chap2.png", surf_2,
       width = the_width, height = the_height, dpi = the_dpi)
Radar Graph
surf_radar_data_2 <- surf %>% 
  filter(climate != "none", chapter %in% c("none", "both","chap_2")) %>% 
  filter(sce != "sce_0", name != "historic run") %>% 
  group_by(timestep, species, code) %>% 
  summarise(kp_nb = mean(kp_nb)) %>% ungroup %>% 
  pivot_wider(names_from = timestep, values_from = kp_nb) %>% 
  mutate(diff = ((`2100`-`2010`)/`2010`)*100) %>% 
  select(-c(`2010`:`2100`)) %>%  
  pivot_wider(names_from = code, values_from = diff) %>% 
  rename(group = species)

surf_radar_2 <- 
  ggradar(surf_radar_data_2, centre.y = -60, legend.position = "right",
          grid.min = -60, grid.max = 5, grid.mid = 0, 
          values.radar = c("-60 %", "0 %", ""), 
          group.colours = group.colours,
          gridline.label.offset = 22, 
          legend.text.size = legend.text.size,
          group.line.width = group.line.width, 
          grid.label.size = grid.label.size,
          background.circle.colour = background.circle.colour, 
          group.point.size  = group.point.size, 
          axis.label.size = axis.label.size-2)
surf_radar_2

ggsave("../thesis/figures/surf_radar_chap2.png", surf_radar_2, 
       width = the_width, height = the_height, dpi = the_dpi)

Maps

list_of_files <-  list.files("../libraries/stsim/monteregie-conncons-scripted.ssim.output/Scenario-37/stsim_OutputSpatialState", full.names = T)
list_of_files <- list_of_files[grep(x=list_of_files, "ts2")]
stratum <- raster("../data/stsim/aggregated/primary_stratum_mont_or_not_or_PA.tif")
names <- paste0(c("historic_compare_agex", "historic_compare_urb"), ".png")

for (class in c(1:2)){
  
  the_stack <- stack(lapply(list_of_files, raster)) == class
  original <- raster("../data/land_use/aggregated/aggregated_lu_buffered_1990_patched.tif") == class
  final <- raster("../data/land_use/aggregated/aggregated_lu_buffered_2010_patched.tif") == class
  
  # Transform, give names
  the_stack[original == 1] <- 0
  final[original == 1] <- 0
  the_stack[stratum==0] <- NA
  final[stratum==0] <- NA
  the_stack <- trim(the_stack)
  final <- trim(final)
  
  the_mean <-  mean(the_stack)
  names(the_mean) <- "mean"
  names(final) <- "observation"
  
  # test <- the_mean
  # test[(which(values(!is.na(final)) & values(is.na(the_mean))))] <- 100
  # plot(test)
  
  # Fill the "other" values
  the_mean[(which(values(!is.na(final)) & values(is.na(the_mean))))] <- 0
  
  stack_to_plot <- stack(the_mean, final) 
  names(stack_to_plot) <- c("Model", "Observations")
  
  mxp <- 1e10000 # 1000
  
  myPal <- RColorBrewer::brewer.pal('OrRd', n=9)
  myTheme <- rasterTheme(region = myPal)
  probs <- levelplot(stack_to_plot$Model, scales=list(draw=FALSE),
                     maxpixels = mxp, 
                     par.settings = myTheme,
                     margin=FALSE, auto.key=FALSE, 
                     colorkey=list(space="bottom"), 
                     main = "A)")
  probs <- as.ggplot(probs) 
  
  stack_to_plot$Observations <- ratify(stack_to_plot$Observations)
  rat <- levels(stack_to_plot$Observations)[[1]]
  rat$Urbanised <- c("Not Observed", "Observed")
  levels(stack_to_plot$Observations) <- rat
  bin <- levelplot(stack_to_plot$Observations, col.regions=c('grey', 'indianred4'), 
                   scales=list(draw=FALSE), 
                   maxpixels = mxp, 
                   colorkey=list(space="bottom", height=0.3), 
                   main = "B)")
  bin <- as.ggplot(bin)
  
  final_plot <- probs + bin
  
  ggsave(plot = final_plot, filename = paste0("../thesis/figures/",names[class]), width = 20, height = 10, dpi = the_dpi)
}

Land use maps

source('../scripts/functions/draw_scenario.R')
for (scenario in c(39, 42, 45, 48, 51)){
  the_plot <- draw_scenario(scenario, mxp = mxp)
  the_name <- paste0("scenario_", as.character(scenario),"_compare.png")
  ggsave(plot = the_plot, 
         filename = paste0("../thesis/figures/",the_name), width = 30, height = 12, dpi = the_dpi)
}

Bar charts

Scenario BAU

final_df_39 <- bar_data %>% 
  filter(scenario == 39)

bar_plot_39 <- final_df_39 %>% 
  #filter(value > 10) %>% 
  mutate(value = as.factor(value), timestep = as.factor(timestep)) %>% 
  ggplot(aes(fill=timestep, y=area, x=new_class)) + 
  geom_bar(position="dodge", stat="identity") +
  #geom_errorbar(aes(ymin=mean_count-sd_count, ymax=mean_count+sd_count), width=.2,
  #               position=position_dodge(.9)) +
  scale_fill_manual(values = bi_colors) +
  theme(legend.position = c(0.85, 0.7)) +
  labs(x = "Land use Class", 
       y = "Mean Area (hectare)", 
       fill = "Year") +
  NULL
bar_plot_39

#ggsave("../thesis/figures/bar_BAU.png", bar_plot_39, width = 15, height = 10)
Scenario Reforestation
# draw_scenario(42)
final_df_42 <- bar_data %>% 
  filter(scenario == 42)

bar_plot_42 <- final_df_42  %>% 
  #filter(value > 10) %>% 
  mutate(value = as.factor(value), timestep = as.factor(timestep)) %>% 
  ggplot(aes(fill=timestep, y=area, x=new_class)) + 
  geom_bar(position="dodge", stat="identity") +
  scale_fill_manual(values = bi_colors) +
  #geom_errorbar(aes(ymin=mean_count-sd_count, ymax=mean_count+sd_count), width=.2,
  #               position=position_dodge(.9)) +
  theme(legend.position = c(0.85, 0.7)) +
  labs(x = "Land use Class", 
       y = "Mean Area (hectare)", 
       fill = "Year") +
  NULL
bar_plot_42

#ggsave("../thesis/figures/bar_Ref.png", bar_plot_42, width = 15, height = 10)
Both scenarios
final_df_39 <- final_df_39 %>% 
  mutate(sce = "BAU") %>% 
  mutate(type = ifelse(value <10, "non_forest", "forest"))
final_df_42 <- final_df_42 %>% 
  mutate(sce = "Reforestation")  %>% 
  mutate(type = ifelse(value <10, "non_forest", "forest"))

final_all <- bind_rows(final_df_39, final_df_42)

forest <- final_all %>% 
  filter(type == "forest") %>% 
  mutate(value = as.factor(value), timestep = as.factor(timestep)) %>% 
  ggplot(aes(fill=timestep, y=area, x=new_class)) + 
  geom_bar(position="dodge", stat="identity") +
  scale_fill_manual(values = bi_colors) +
  theme(legend.position = "none") +
  facet_grid(~sce) +
  labs(x = "Land use Class", 
       y = "Mean Area (hectare)", 
       fill = "Year") +
  NULL

non_forest <- final_all %>% 
  filter(type != "forest") %>% 
  filter(!(new_class %in% c("Water", "Wetlands", "Roads"))) %>% 
  mutate(value = as.factor(value), timestep = as.factor(timestep)) %>% 
  ggplot(aes(fill=timestep, y=area, x=new_class)) + 
  geom_bar(position="dodge", stat="identity") +
  scale_fill_manual(values = bi_colors) +
  theme(legend.position = "right") +
  facet_grid(~sce) +
  labs(x = "Land use Class", 
       y = "Mean Area (hectare)", 
       fill = "Year") +
  NULL

assembled <- wrap_plots(forest, non_forest, ncol=1, nrow=2)
assembled

ggsave("../thesis/figures/bar_both.png", assembled,
       width = the_width, height = 20, dpi = the_dpi)
Other calculations

Compute auc / precision on predictions.